Iterative calibration method for a direct neural interface using a markov mixture of experts with multivariate regression

ABSTRACT

This invention relates to a method of calibrating a direct neural interface with continuous coding. The observation variable is modelled by an HMM model and the control variable is estimated by means of a Markov mixture of experts, each expert being associated with a state of the model. 
     During each calibration phase, the predictive model of each of the experts is trained on a sub-sequence of observation instants corresponding to the state with which it is associated, using an REW-NPLS (Recursive Exponentially Weighted N-way Partial Least Squares) regression model. 
     A second predictive model giving the probability of occupancy of each state of the HMM model is also trained during each calibration phase using an REW-NPLS regression method. This second predictive model is used to calculate Markov mixture coefficients during a later operational prediction phase.

TECHNICAL FIELD

This invention relates to the field of direct neural interfaces alsoknown as BCI (Brain Computer Interface) or BMI (Brain MachineInterface). Its applications are particularly direct neural control of amachine such as an exoskeleton or a computer.

STATE OF PRIOR ART

Direct neural interfaces use electrophysiological signals emitted by thecerebral cortex to generate a control signal. Much research has beenmade on these neural interfaces, particularly with the purpose ofrestoring a motor function to a paraplegic or tetraplegic patient usinga prosthesis or a motor-driven orthosis.

Neural interfaces can have an invasive or a non-invasive nature.Invasive neural interfaces use intracortical electrodes (in other wordsimplanted in the cortex) or cortical electrodes (positioned on thesurface of the cortex), in the latter case collectingelectrocorticographic (ECoG) signals. Non-invasive neural interfaces useelectrodes placed on the scalp to collect electroencephalographic (EEG)signals. Other types of sensors have also been envisaged such asmagnetic sensors measuring magnetic fields induced by the electricalactivity of brain neurons. The term used in this case ismagnetoencephalographic (MEG) signals.

Direct neural interfaces advantageously use ECoG type signals that havethe advantage of providing a good compromise between biocompatibility(matrix of electrodes implanted on the surface of the cortex) andquality of the collected signals.

The ECoG signals thus measured must be processed so as to estimate themovement trajectory required by the patient and deduce control signalsfor the computer or the machine from these signals. For example, whenthe objective is to control an exoskeleton, the BCI interface estimatesthe required movement trajectory from the measured electrophysiologicalsignals and uses them to deduce the control signals that will enable toexoskeleton to reproduce the trajectory in question. Similarly, when theobjective is to control a computer, the BCI interface may for exampleestimate the required trajectory of a pointer or a cursor starting fromelectrophysiological signals and use them to deduce cursor/pointercontrol signals.

The estimate of the trajectory, and more precisely of the kinematicparameters (position, speed and acceleration), is also known as neuraldecoding in the literature. In particular, neural decoding can be usedto control a movement (of a prosthesis or a cursor) from ECoG signals.

When the ECoG signals are acquired continuously, one of the maindifficulties in decoding lies in the asynchronous nature of the control,in other words in the discrimination of phases during which the patientactually controls a movement (active period) and phases in which he doesnot control it (idle periods).

To circumvent this difficulty, direct neural interfaces calledsynchronous interfaces can only control a movement during well-definedtime windows (for example periodically succeeding time intervals)signalled to the patient by an external index. The patient can then onlycontrol the movement during these time windows, but this proves to beunworkable in most practical applications.

More recently, direct neural interfaces with continuous decoding havebeen disclosed in the literature. The paper by M. Velliste et al.entitled “Motor cortical correlates of arm resting in the context of areaching task and implications for prosthetic control” published in TheJournal of Neuroscience, Apr. 23, 2014, 34 (17), pp. 6011-6022,describes in particular a direct neural interface in which idle stateperiods are detected by LDA (Linear Discriminant Analysis) making use ofthe action potentials emission frequency (neuron firing rate). Kinematicparameters are estimated during active periods by means of a statetransition model, the states being predicted by a Laplacian of Gaussianfilter. However, the results have been obtained for matrices ofmicro-electrodes and could not be reproduced for classical corticalelectrodes. Furthermore, switchings between idle periods and activeperiods result in discontinuities in decoding the movement and thereforesudden changes along the trajectory.

Another approach has been proposed in the paper by L. Srinivasam et al.entitled “General-purpose filter design for neural prosthetic devices”published in J. Neurophysiol. Vol 98, pp. 2456-2475, 2007. This paperdescribes a direct neural interface with continuous decoding in whichthe sequence of active states and idle states is generated by a 1^(st)order Markov model. A Switching Kalman Filter or SKF for which theobservation matrices and the transition matrices depend on the hiddenswitching variable (active state/idle state) can be used to estimate thekinematic parameters of the movement. This type of direct neuralinterface has been used to control a wheelchair making use of EEGsignals based on simulation data but has not been tested for ECoGsignals. Furthermore, detection of the state (active state/idle state)is sometimes incorrect (high rate of false positives and falsenegatives).

Finally, the above-mentioned direct neural interfaces have beendeveloped to decode the movement of a single member of a patient.Therefore they are not adapted to the control of an exoskeleton withseveral members, particularly for a tetraplegic, paraplegic orhemiplegic patient, which significantly reduces their application field.

In the application published under number FR-A-3053495, it is proposedto decode the movement of several members by using a Markov mixture ofexperts or MME. This decoding is based on a Hidden Markov Model or HMM,an estimator (or expert) being associated with each hidden state of themodel.

FIG. 1 schematically represents a BCI interface using a continuousmovement decoding by a Markov mixture of experts.

This interface, 100, makes use firstly of a hidden state machine, 140,that can take K possible states, and secondly, a plurality K ofestimators, 120, each estimator (or expert) being associated with ahidden state.

Estimates made by the different experts are combined in a combinationmodule (gating network), 130, to give an estimate, ŷ(t) of the variableto be explained, y(t), in this case the kinematic movement parameters,starting from the observation x(t) representing the characteristics ofneural signals at time t, captured by means of electrodes 105. Thesequence of observations x[1:t]={x(1), x(2), . . . , x(t)} is modelledby a subjacent first order Markov process with continuous emission.

If the input variable x(t) has dimension M and the response y(t) hasdimension N, and if the predictive models of the different expertsE_(k), are chosen to be linear, in other words E_(k)(x(t))=β_(k)x(t) inwhich β_(k) is a matrix with size N×M, the estimate by Markov mixture ofexperts is expressed as follows:

$\begin{matrix}{{\overset{\hat{}}{y}(t)} = {\sum\limits_{k = 1}^{K}{{g_{k}(t)}\beta_{k}{x(t)}}}} & (1)\end{matrix}$

in which g_(k)(t) are weighting (or mixing) coefficients satisfying

${\sum\limits_{k = 1}^{K}{g_{k}(t)}} = 1.$

As indicated in the above-mentioned application, the parameters of theHMM model, in other words the probabilities of transition between hiddenstates, and the matrices β_(k) can be estimated during a calibrationphase, making use of a supervised estimating method. In this case, theHMM machine and the K experts are trained independently of each other.In particular, each expert E_(k) can be trained on a sub-sequence ofobservations corresponding to the state k, the matrix β_(k) of thepredictive model then being estimated using a linear PLS (Partial LeastSquares) regression starting from this sub-sequence.

More generally, the linear relation (1) can be extended to amulti-linear relation. The input variable can then be represented intensor form X∈□^(I) ¹ ^(× . . . ×I) ^(n) n in which n is the order ofthe input tensor and I_(i) is the dimension of mode i.

Input tensor modes can for example be time (number of samples in time),frequency (number of spectral analysis bands), space (number ofelectrodes). Similarly, the response may be represented in tensor formY∈□^(J) ¹ ^(× . . . ×J) ^(m) in which in is the order of the outputtensor. Output tensor modes can correspond to different effectors, forexample acting on different joints of an exoskeleton.

Tensors of predictive models of the different experts can then beestimated by a multivariate PLS or NPLS (N-way PLS) method as describedin the application published under number FR-A-3046471.

Regardless of which predictive models (PLS or NPLS) of the differentexperts are used, it has been observed that the validity duration ofthese models is relatively limited due to the lack of stationarity ofsignals in the human brain. Consequently, they have to be updatedregularly by performing new calibration phases. However, a very largequantity of data has to be processed for each new calibration phase,combining data from previous calibrations and from the new calibration,for the different experts. The duration of the update is then often suchthat functioning of the direct neural interface has to be interrupted(off-line calibration). An iterative method of calibrating a directneural interface based on an NPLS predictive model has thus beenproposed in application FR-A-3061318. This iterative calibration methodwas also described in the paper written by A. Eliseyev et al. entitled“Recursive exponentially weighted N-way Partial Least Squares regressionwith recursive validation of hyper-parameters in Brain ComputerInterface applications” published in Nature, Scientific Reports, volume7, article number 16281, published on 24 Nov. 2017.

However, this iterative calculation method does not apply to a directneural interface using a Markov mixture of experts as described inapplication FR-A-3046471.

One purpose of this invention is consequently to disclose a new methodof calibrating a direct neural interface using a Markov mixture ofexperts, without needing to process a considerable quantity of dataduring each calibration phase, thus making operation of the directneural interface possible in real time, without interruption.

PRESENTATION OF THE INVENTION

This invention is defined by a method of calibrating a direct neuralinterface according to the description in claim 1 in the set of claimsgiven below. Advantageous embodiments are specified in the dependentclaims.

BRIEF DESCRIPTION OF THE DRAWINGS

Other characteristics and advantages of the invention will become clearafter reading a preferred embodiment of the invention, described withreference to the appended figures among which:

FIG. 1, described above, schematically represents a direct neuralinterface using a Markov mixture of experts, known in prior art;

FIG. 2 schematically represents a direct neural interface using a Markovmixture of experts according to one embodiment of the invention;

FIG. 3 is a flowchart representing a method of iteratively calibrating adirect neural interface using a Markov mixture of experts according toone embodiment of the invention.

DETAILED PRESENTATION OF PARTICULAR EMBODIMENTS

In the following, we will consider a direct neural interface withcontinuous decoding using a Markov mixture of experts as described inthe introductory part.

The electrophysiological signals output from the different electrodesare sampled and assembled by data blocks, each block corresponding to asliding observation window with width ΔT. For example, theelectrophysiological signals are sampled at a frequency of about 590kHz, and the size of data blocks is 590 samples (ΔT=1 s), and the offsetδT from one observation window to the next is 59 samples (δT=100 ms).Each observation window is defined by an observation time (epoch) atwhich the window in question starts.

The electrophysiological signals are then preprocessed. In particular,this preprocessing may comprise elimination of the average taken on allelectrodes, and a time-frequency analysis is then made on eachobservation window.

The time-frequency analysis can use a decomposition into wavelets,particularly Morlet wavelets.

A frequency smoothing or decimation can also be made on these results ofthe time-frequency analysis.

Thus, each observation window, or observation instant t, is associatedwith an order 3 tensor of observation data, resulting in the generationof an order 4 input tensor: the first mode corresponds to successiveobservation windows, the second mode corresponds to space, in otherwords the sensors, the third mode corresponds to time within anobservation window, in other words the positions of the wavelets, andthe fourth mode corresponds to the frequency, in other words to thenumber of frequency bands used for the decomposition into wavelets on anobservation window.

More generally, the input tensor (or the observation tensor) will beorder n+1, the first mode always being the mode relative to observationtimes (epochs). The input tensor (or observation tensor) is denoted Xand its dimension is N×I₁× . . . ×I_(n).

Similarly, the trajectory of the imagined movement, observed orperformed, is described by an output tensor (or command tensor), orderm+1, denoted Y, with dimension N×J₁× . . . ×J_(m), of which the firstmode corresponds to successive times at which the commands will beapplied (in general, this first mode also corresponds to the observationwindows), the other modes corresponding to commands on the differenteffectors or different degrees of freedom of a multi-axis robot.

More precisely, the output tensor supplies N consecutive command datablocks, each of the blocks being used to generate command signals forthe different effectors or degrees of freedom. Thus, it will beunderstood that the dimension of each data block can depend on theenvisaged usage case and particularly on the number of degrees offreedom of the effector.

In the following, X _(t) will be used to denote the observation tensorat time t. Consequently, this tensor is order n and its dimension is I₁×. . . ×I_(n). It takes its values in a space X⊂□^(I) ¹ ^(× . . . ×I)^(n) in which □ is the set of reals. Similarly, Y _(t) will be used todenote the command tensor at time t. Consequently, this output tensor isorder m and dimension J₁× . . . ×J_(m). It takes its values in a spaceY⊂□^(J) ¹ ^(× . . . ×J) ^(m) .

It will be assumed that the space X is formed from the union of aplurality K of regions, not necessarily discontiguous. In other words,

$\underset{¯}{X} = {\bigcup\limits_{k = 1}^{K}X_{k}}$

in which X_(k), k=1, . . . , K are the regions in question.

The direct neural interface is based on a mixture of experts (ME), eachexpert E_(k) operating on the elementary region X_(k) of the inputcharacteristics space and being capable of predicting the command tensorY _(t) at time t starting from the observation tensor, X _(t), when thistensor belongs to X_(k). In other words, each region X_(k) has acorresponding command tensor prediction model Y _(t) starting from theobservation tensor X _(t). An expert E_(k) can thus be considered as amulti-linear application of X_(k) in Y.

It will also be assumed that each expert E_(k) is associated with ahidden state k of a first order Markov Model (HMM). The differentexperts are combined by means of combination coefficients dependent onthe hidden state at the time of the prediction. Definitively, startingfrom an input (or observation) tensor X _(t), the neural interfaceestimates the output (or command) tensor Y _(t), using:

$\begin{matrix}{{\bullet {\underset{\_}{Y}}_{t}} = {\sum\limits_{k = 1}^{K}{\gamma_{k}^{t}\left( {{{\underset{\_}{\beta}}_{k}{\underset{\_}{X}}_{t}} + {\underset{\_}{\delta}}_{k}} \right)}}} & (2)\end{matrix}$

in which β _(k), δ _(k) are a prediction coefficients tensor and anexpert prediction bias tensor respectively, E_(k) and δ_(k) ^(t) is theweighting coefficient (also called the gating coefficient) for thisexpert and at the time t considered. The weighting coefficient γ_(k)^(t) is simply the conditional probability that the HMM model is in thestate k knowing the previous input tensors X _(1:t)=X ₁, X ₂, . . . , X_(t).

The set of coefficients and prediction biases, collectively calledprediction parameters of the different experts, is designated byθ_(e)={(β _(k),δ_(k))|k=1, . . . , K}. The set of parameters for thecombination of the different experts (including the parameters of thesubjacent HMM model) is designated by θ_(g)={A,{d_(k)|k=1, . . . , K},π} in which A is the transition matrix between states, with size K×K, inother words the matrix for which the elements a^(ij) (independent oftime by assumption of the HMM model) are defined bya^(ij)=p(z_(t)=j|z_(t−1)=i) in which z, represents the state of themodel at time t and z_(t−1) represents the state of the model at thepreceding time; {d_(k)|k=1, . . . , K} the parameters used to determinethe conditional probability of observing the input tensor X _(t) knowingthe state z_(t)=k, and π is a vector with size K giving probabilities ofoccupancy of the different states at the initial time, in other wordsπ_(i)=p(z_(t)=i)_(t=0).

FIG. 2 schematically represents a direct neural interface using a Markovmixture of experts according to one embodiment of the invention.

The electrophysiological signals output from the different electrodes205 are processed in the processing module 210. This processing moduleperforms sampling, optionally eliminates the average taken on allelectrodes, and a time-frequency analysis is then made on eachobservation window. The time-frequency analysis can possibly be followedby frequency smoothing or decimation, as indicated above. Consequently,the processing module 210 outputs an input tensor X _(t) at eachobservation window, characterised by an observation time t.

The direct neural interface also makes use of a hidden state machinebased on an HMM model, 240, that can be in K possible states, and aplurality of experts, 220, each expert being associated with a hiddenstate.

Each expert E_(k) makes an estimate Y _(t) ^(k) of the output vector Y_(t) and supplies it to a combination module, 230, that mixes theestimates according to expression (2). The combination coefficients(gating coefficients) are determined in the estimating module, 250, byγ_(k) ^(t)=p(z_(t)=k|X _(1:t)), k=1, . . . , K, in which z_(t)represents the state of the model at time t.

The estimate

_(t) gives the different kinematic movement parameters such as theposition and/or the speed and/or the acceleration along the requiredtrajectory.

The neural interface is calibrated firstly during an initialisationphase and then, during subsequent calibration phases, at given timeintervals (on-line calibration).

These subsequent calibration phases are necessary due to the lack ofstationarity of signals generated in the human brain. They can be doneat regular intervals along the trajectory or successive trajectories tobe executed.

A calibration phase u uses a plurality of input data blockscorresponding to a plurality ΔL of successive times output from theobservation times sequence. X _(u) will denote the tensor with dimensionΔL×I₁× . . . ×I_(n) containing this plurality of input blocks. It alsomakes use of a plurality of output blocks giving kinematic setparameters, taken at the same times. Y _(u) will denote the set tensorwith dimension ΔL×J₁× . . . ×J_(m) containing this plurality of inputblocks.

The switches 261, 262 are in the closed position during a calibrationphase. The calibration module 270 receives the input tensor X _(u), theset tensor Y _(u), and the matrix Z_(u) with size K×ΔL, of which each ofthe columns represents the state of the machine at the time considered.More precisely, if at time

the machine is in state z

=k, all elements in the

^(th) column of the matrix Z_(u) are null except for the element on rowk that is equal to 1.

The calibration module uses X _(u), Y _(u), Z _(u) to calculate theupdated parameter sets θ_(e) and θ_(g).

At the end of this calibration phase, the switches 261, 262 are open andthe updated parameter sets θ_(e) and θ_(g) are output to a plurality ofexperts, 220, and to the machine 240, respectively.

Innovatively, the predictive model of each expert E_(k) is trained usingan REW-NPLS (Recursive Exponentially Weighted N-way Partial LeastSquares) regression with a forget factor λ^(k) (0<λ^(k)<1) oncalibration data blocks, X _(u), Y _(u), corresponding to the hiddenstate z_(u)=k, denoted X _(u) ^(k) and Y _(u) ^(k). More precisely, foreach expert E_(k) and for each calibration phase u, the data blocks X_(u) ^(k) and Y _(u) ^(k) are determined so that it can be trained in asupervised manner.

This is done as described in application FR-A-3061318, namely:

In a first step, the tensors X _(u) ^(k) and Y _(u) ^(k) are normalisedstarting from a number of observations said to be effective N_(u)^(k)=λ^(k)N_(u−1) ^(k)+N_(u) ^(k) in which N_(u) ^(k) is the number ofdata blocks corresponding to the hidden state z_(u)=k during thecalibration phase u. The weighted sum s(X _(u) ^(k)) and the weightedquadratic sum sq(X _(u) ^(k)) are calculated:

$\begin{matrix}{{s\left( {\underset{\_}{X}}_{u}^{k} \right)} = {{\lambda^{k}{\sum\limits_{ = 1}^{N_{u - 1}^{k}}x_{{{u - 1};},i_{1},\; {\ldots \mspace{11mu} i_{n}}}^{k}}} + {\sum\limits_{ = 1}^{N_{u}^{k}}x_{{u;},i_{1},\; {\ldots \mspace{11mu} i_{n}}}^{k}}}} & \left( {3\text{-}1} \right) \\{{{sq}\left( {\underset{\_}{X}}_{u}^{\;_{k}} \right)} = {{\lambda^{k}{\sum\limits_{ = 1}^{N_{u - 1}}\left( x_{{{u - 1};},i_{1},\; {\ldots \mspace{11mu} i_{n}}}^{k} \right)^{2}}} + {\sum\limits_{ = 1}^{N_{u}^{k}}\left( x_{{u;},i_{1},\; {\ldots \mspace{11mu} i_{n}}}^{k} \right)^{2}}}} & \left( {3\text{-}2} \right)\end{matrix}$

in which

and

represent the elements of tensors X _(u−1) ^(k) and X _(u) ^(k).respectively. The average value

${\mu \left( {\underset{\_}{X}}_{u}^{k} \right)} = \frac{s\left( {\underset{\_}{X}}_{u}^{k} \right)}{N_{u}^{k}}$

and the standard deviation

${\sigma \left( {\underset{¯}{X}}_{u}^{k} \right)} = \sqrt{\frac{{s{q\left( {\underset{¯}{X}}_{u}^{k} \right)}} - \left( {\mu \left( {\underset{¯}{X}}_{u}^{k} \right)} \right)^{2}}{N_{u}^{k}}}$

are then deduced. The elements of the centred and normalised inputtensor,

, are:

$\begin{matrix}{{\overset{\sim}{x}}_{{u;},i_{1},\; {\ldots \mspace{11mu} i_{n}}}^{k} = \frac{x_{{u;},i_{1},\; {\ldots \mspace{11mu} i_{n}}}^{k} - {\mu \left( {\underset{\_}{X}}_{u}^{k} \right)}}{\sigma \left( {\underset{\_}{X}}_{u}^{k} \right)}} & (4)\end{matrix}$

Similarly, the weighted sum s(Y _(u) ^(k)) and the weighted quadraticsum sq(Y _(u) ^(k)) are calculated:

$\begin{matrix}{{s\left( {\underset{¯}{Y}}_{u}^{k} \right)} = {{\lambda^{k}{\sum\limits_{ = 1}^{N_{u - 1}^{k}}y_{{{u - 1};},j_{1},\; {\ldots \mspace{11mu} j_{m}}}^{k}}} + {\sum\limits_{ = 1}^{N_{u}^{k}}y_{{u;},j_{1},\; {\ldots \mspace{11mu} j_{m}}}^{k}}}} & \left( {5\text{-}1} \right) \\{{{sq}\left( {\underset{¯}{Y}}_{u}^{k} \right)} = {{\lambda^{k}{\sum\limits_{ = 1}^{N_{u - 1}^{k}}\left( y_{{{u - 1};},j_{1},\; {\ldots \mspace{11mu} j_{m}}}^{k} \right)^{2}}} + {\sum\limits_{ = 1}^{N_{u}^{k}}\left( y_{{u;},j_{1},\; {\ldots \mspace{11mu} j_{m}}}^{k} \right)^{2}}}} & \left( {5\text{-}2} \right)\end{matrix}$

The average value

${\mu \left( {\underset{\_}{Y}}_{u}^{k} \right)} = \frac{s\left( {\underset{¯}{Y}}_{u}^{k} \right)}{N_{u}^{k}}$

and the standard deviation

${\sigma \left( {\underset{¯}{Y}}_{u}^{k} \right)} = \sqrt{\frac{{{sq}\left( {\underset{¯}{Y}}_{u}^{k} \right)} - \left( {\mu \left( {\underset{¯}{Y}}_{u}^{k} \right)} \right)^{2}}{N_{u}^{k}}}$

are then deduced, followed by the elements of the centred and normalisedoutput tensor,

:

$\begin{matrix}{{\overset{\sim}{y}}_{{u;},j_{1},\; {\ldots \mspace{11mu} j_{m}}}^{k} = \frac{y_{{u;},j_{1},\; {\ldots \mspace{11mu} j_{m}}}^{k} - {\mu \left( {\underset{¯}{Y}}_{u}^{k} \right)}}{\sigma \left( {\underset{¯}{Y}}_{u}^{k} \right)}} & (6)\end{matrix}$

The covariance tensor and the cross-covariance tensor are defined duringthe calibration phase u, starting from the centred and normalisedtensors

and

:

$\begin{matrix}{{{cov}(,)} = { \times_{1}}} & \left( {7\text{-}1} \right) \\{{{cov}(,)} = { \times_{1}}} & \left( {7\text{-}2} \right)\end{matrix}$

in which x₁ designates the tensor product according to the first mode(temporal mode with index

) as described in the above-mentioned application FR-A-3061318.

These covariance tensors are modified taking account of the covariancetensors of the previous calibration phase, by weighting them with theforget factor λ^(k):

cov(

,

)=λ^(k) cov(

,

)+

x ₁

  (8-1)

cov(

,

)=λ^(k) cov(

,

)+

x ₁

  (8-2)

This calculation takes account of the data from a previous calibrationphase, with a forget factor, to update the covariance tensors. Theforget factors relative to the different experts can be chosen to beidentical λ^(k)=λ, k=1, . . . , K, in which λ is then the common forgetfactor. Alternatively, they can be chosen to be distinct so as to offermore prediction flexibility to the different experts.

Starting from the covariance tensors cov(

,

) and cov(β,

) working iteratively on rank f^(k)=1, . . . , F in the latent variablesspace (in which F is a maximum predetermined number of dimensions in thelatent variables space), we obtain a set of projection vectors {w_(u,1)^(k,f) ^(k) , . . . , w_(u,n) ^(k,f) ^(k) }_(f) _(k) ₌₁ ^(F) withdimensions I₁, . . . , I_(n) respectively, a set of predictioncoefficient tensors, {β _(u) ^(k,f) ^(k) }_(f) _(k) ₌₁ ^(F), and aprediction bias sensor, {δ _(u) ^(k,f) ^(k) }_(f) _(k) ₌₁ ^(F), usingthe REW-NPLS method described in the above-mentioned applicationFR-A-3061318.

The optimum rank, f_(opt) ^(k), for each expert E_(k), during acalibration phase u is determined by calculating the prediction errors Y_(u) ^(k) starting from X _(u) ^(k) using the prediction coefficienttensors, {β _(u−) ^(k,f) ^(k) }_(f) _(k) ₌₁ ^(F), and the predictionbias tensors, {δ _(u−) ^(k,f) ^(k) }_(f) _(k) ₌₁ ^(F) obtained duringthe previous calibration phase:

err_(u) ^(k,f) ^(k) =λ^(k) err_(u−1) ^(k,f) ^(k) +∥Y _(u) ^(k) −Ŷ ^(k,f)^(k) ∥  (9)

in which Ŷ ^(k,f) ^(k) is the estimate of Y _(u) ^(k) obtained from theprediction β _(u−1) ^(k,f) ^(k) and prediction bias coefficients δ_(u−1) ^(k,f) ^(k) from the previous calibration phase.

The rank f_(opt) ^(k) leading to the minimum prediction error is thenselected:

$\begin{matrix}{f_{opt}^{k} = {\underset{{f^{k} = 1},\; \ldots \;,F}{\arg \; \min}\left( {err}_{u}^{k,f^{k}} \right)}} & (10)\end{matrix}$

and for the calibration phase, the prediction coefficients tensor andthe prediction bias tensor corresponding to this optimal rank areselected, in other words the update is made at the end of thecalibration phase u:

β _(k)=β _(u) ^(k,f) ^(opt) ^(k)   (11-1)

δ _(k)=δ _(u) ^(k,f) ^(opt) ^(k)   (11-2)

for each of the experts E_(k), k=1, . . . , K, independently. The resultobtained is thus an estimate of the set of prediction parameters θ_(e)making use of the REW-NPLS regression method applied to each of theexperts.

Innovatively, the combination coefficients γ_(k) ^(t), k=1, . . . , Kare also estimated using an REW-NPLS regression method adapted todiscrete decoding, as described later.

During a calibration phase, the elements a^(ij) of the transition matrixA are firstly updated as follows:

$\begin{matrix}{a_{u}^{ij} = {{\lambda a_{u - 1}^{ij}} + \frac{v_{u}^{ij}}{K}}} & (12)\end{matrix}$

in which a_(u) ^(ij) (or a_(u−1) ^(ij)) are elements of the transitionmatrix, estimated during the calibration phase u (or u−1) and is thenumber of transitions from state i to state j during the calibrationphase u. This number of transitions is obtained starting from the matrixZ_(u). The lines in the transition matrix are then normalised such thesum of the probabilities of the transition from an initial state isequal to 1.

The transition matrix A thus updated is supplied by the transitionmatrix calculation module, 260, to the state estimating module, 240, ofthe HMM model.

z_(l) defines a vector with size K for which the elements representwhether or not the machine is in state k=1, . . . , K at time t. Thus,if the machine is in state k at time t, only the kth element of z_(t)will be equal to 1 and the other elements will be null. z_(t) can beconsidered like a process with discrete values that vary in timeaccording to a multilinear predictive model, the input variable of whichis the observation tensor X _(t). This predictive model is trained undersupervision during each calibration phase u, using the same principle asthe predictive model of each expert. In other words:

{circumflex over (z)} _(t) =B X _(t) +b   (13)

in which {circumflex over (z)}_(t) expresses the probabilities{circumflex over (z)}_(k,t) that the machine is in state k at time t, Bis a prediction coefficients tensor with size K×I₁× . . . ×I_(n) and bis a bias vector with size K.

The tensor B and the vector b are updated during each calibration phaseu by means of an REW-NPLS regression with a forget factor i λ, (0<λ<1)based on observations X _(u) and Z_(u) (matrix with size K×ΔL). Thisforget factor may be identical to the common forget factor, when such acommon forget factor is used to estimate prediction parameters for thedifferent experts.

Working iteratively on rank f=1, . . . , F of the latent variables spacestarting from the covariance tensor cov(

,

) in which

is a centred and normalised version of X _(u), and the cross-covariancetensor cov(

, Z_(u)). These covariance and cross-covariance tensors are modified by:

cov(

,

)=λ·cov(

,

)+

x ₁

  (14-1)

cov(

,Z_(u))=λ·cov(

,Z_(u))+

x ₁ Z _(u)   (14-2)

The REW-NPLS method provides a set of projection vectors {w_(u,1), . . ., w_(u,n)}_(f=1) ^(F), with dimensions I₁, . . . , I_(n) respectively, aset of prediction coefficient tensors, {β _(u) ^(f)}_(f=1) ^(F), and aset of prediction bias vectors, {b_(u) ^(f)}_(f=1) ^(F). Rank f_(opt)leading to the minimum prediction error is selected, then the tensor Band the vector b are updated by:

B=B _(u) ^(f) ^(opt)   (15-1)

b=b_(u) ^(f) _(opt)   (15-2)

The tensor B and the vector b are supplied to the mixture coefficientsestimating module, 250.

Starting from elements {circumflex over (z)}_(k,t) of the predictedvector {circumflex over (z)}_(t), the state estimating module cancalculate the conditional probabilities using the softmax function:

$\begin{matrix}{{p\left( {z_{t} = \left. k \middle| {\underset{¯}{X}}_{t} \right.} \right)} = \frac{\exp \left( {\overset{\hat{}}{z}}_{k,t} \right)}{\sum\limits_{i = 1}^{K}{\exp \left( {\overset{\hat{}}{z}}_{i,t} \right)}}} & (16)\end{matrix}$

This expression gives the probability that the machine is in state k attime t, knowing the input data at this time.

The mixing coefficients estimating module, 250, obtains thesecoefficients from Bayes rule:

$\begin{matrix}{\gamma_{k}^{t} = {{p\left( {{z_{t} = \left. k \middle| {\underset{\_}{X}}_{1:t} \right.},} \right)} = {\frac{p\left( {{z_{t} = k},{\underset{\_}{X}}_{1:t}} \right)}{p\left( {\underset{\_}{X}}_{1:t} \right)} = \frac{p\left( {{z_{t} = k},{\underset{\_}{X}}_{1:t}} \right)}{\sum\limits_{i = 1}^{K}{p\left( {{z_{t} = i},{\underset{\_}{X}}_{1:t}} \right)}}}}} & (17)\end{matrix}$

in which p(z_(t)=i, X _(1:t)) are obtained using the forward algorithm,namely:

$\begin{matrix}{{p\left( {{z_{t} = i},{\underset{\_}{X}}_{1:t}} \right)} = {{p\left( {\left. {\underset{\_}{X}}_{t} \middle| z_{t} \right. = i} \right)} \cdot {\sum\limits_{\;^{j = 1}}^{K}{a^{ij}\gamma_{k}^{t - 1}}}}} & (18)\end{matrix}$

The mixture coefficients can be obtained by recurrence by combiningexpressions (17) and (18):

$\begin{matrix}{\gamma_{k}^{t} = \frac{{p\left( {\left. {\underset{\_}{X}}_{t} \middle| z_{t} \right. = k} \right)} \cdot {\sum\limits_{\;^{j = 1}}^{K}{a^{kj}\gamma_{k}^{t - 1}}}}{\sum\limits_{i = 1}^{K}\left( {{p\left( {\left. {\underset{\_}{X}}_{t} \middle| z_{t} \right. = i} \right)} \cdot {\sum\limits_{j = 1}^{K}{a^{ij}\gamma_{i}^{t - 1}}}} \right)}} & (19)\end{matrix}$

The conditional emission probabilities, p(X _(t)|z_(t)=i) can beobtained using Bayes rule starting from a posteriori conditionalprobabilities p(z_(t)=i|X _(t)) and a priori probabilities p(z_(l)=i|X_(t)):

$\begin{matrix}{{p\left( {\left. {\underset{\_}{X}}_{t} \middle| z_{t} \right. = i} \right)} = {\frac{{p\left( {z_{t} = \left. i \middle| {\underset{\_}{X}}_{t} \right.} \right)}{p\left( {\underset{\_}{X}}_{t} \right)}}{p\left( {z_{t} = i} \right)} \propto \frac{p\left( {z_{t} = \left. i \middle| {\underset{\_}{X}}_{t} \right.} \right)}{p\left( {z_{t} = i} \right)}}} & (20)\end{matrix}$

the term p(X _(t)) being a multiplication coefficient common to allstates.

The probabilities of occupancy of the different states p(z_(t)=i) attime t are calculated from the initial probabilities of occupancy π andthe transition matrix A.

FIG. 3 is a flowchart representing a method of iteratively calibrating adirect neural interface using a Markov mixture of experts according toone embodiment of the invention;

Said calibration method works during calibration phases planned atdetermined instants, for example at regular time intervals throughoutthe trajectory. Each calibration phase u uses a plurality of input datablocks corresponding to a plurality ΔL of successive times in theobservation times sequence. It also makes use of a plurality of outputblocks giving kinematic set parameters at the same times.

The plurality of input data blocks for the calibration phase u isrepresented by the input tensor X _(u) and the plurality of outputblocks during this same calibration phase is represented by the outputtensor Y _(u).

It is also assumed that the states of the HMM machine are known duringthe calibration phase. The different states of the HMM machine canrelate to different elements to be controlled to make the trajectory,for example different members of an exoskeleton. Furthermore, for eachof these elements, different states can be considered for example suchas an active state when the patient controls this element, an idle statewhen he does not control it, and possibly a preparation stateimmediately preceding the active state and immediately following theidle state, during which the characteristics of the neural signals arenot the same as in the idle state nor the active state for this element.Thus, when the direct neural interface must control P elements to makethe trajectory, the HMM machine will comprise 2^(P) or even 3^(P) states(if preparatory states are envisaged). The states of the HMM machineduring the calibration phase can be represented by a binary matrix Z_(u)with size K×ΔL, all values in each column in Z_(u) being 0 except forone element equal to 1 indicating the state of the machine at thecorresponding time.

The input tensor, X _(u), the output tensor, Y _(u), the states matrix,Z_(u), are given in 310 for the calibration phase u.

In step 320, the matrix Z_(u) is used as a starting point to determinethe tensors X _(u) ^(k) and Y _(u) ^(k) formed by the input data blocksand the output data blocks related to the state k=1, . . . , K,respectively. These tensors are extracted from the input tensor, X _(u),and the output tensor, Y _(u).

In step 330, for each expert E_(k), k=1, . . . , K the centred andnormalised tensors

and

are calculated and are used to deduce the covariance tensor

, and the cross-covariance tensor of

and

. These tensors are modified by adding the covariance tensor of

⁻¹ and the cross-covariance tensor respectively, obtained in theprevious calibration phase, weighted by a forget factor λ^(k). Thesetensors thus modified will be used as covariance and cross-covariancetensors during the next calibration phase.

In step 340, the prediction parameter tensors β _(k), δ _(k) of eachexpert E_(k) are updated using an REW-NPLS regression, starting from thecovariance and cross-covariance tensors modified in the previous step.

After step 340, we have an updated version of the set θ_(e) of expertprediction parameters K.

In step 350, the elements of the transition matrix A are advantageouslyupdated starting from the number of transitions between successivestates observed during the calibration phase u. The number oftransitions between successive states is obtained starting from thematrix Z_(u). The matrix A is then modified by adding to it the matrix Aobtained during the previous iteration, weighted by a forget factor λ.

In step 360, the centred and normalised tensor

is added, followed by the covariance tensor of

and the cross-covariance tensor of

and Z_(u). These tensors are modified by adding the covariance tensor of

and the cross-covariance tensor

and Z_(u−1), respectively, obtained in the previous calibration,weighted by the forget factor λ. These tensors thus modified will beused as covariance and cross-covariance tensors during the nextcalibration phase.

In step 370, a multi-linear predictive model is trained giving a statevector of the HMM machine as a function of the input tensor, {circumflexover (z)}_(t)=B X _(t)+b, the components of the state vector {circumflexover (z)}_(t) providing probabilities that the machine is in thedifferent states k=1, . . . , K respectively at time t. More precisely,the prediction parameter tensor B, and the bias vector b are updatedusing an REW-NPLS regression, starting from the covariance andcross-covariance tensors modified in the previous step.

The components of {circumflex over (z)}_(t) calculated during anoperational phase are then used at each observation time t to obtainmixture coefficients γ_(k) ^(t)=p(z_(t)=k|X _(1:t)) of the differentexperts.

1. A method of calibrating a direct neural interface that will receive aplurality of electrophysiological signals acquired by means of aplurality of sensors, during a plurality of observation windowsassociated with observation times, and to provide command signalsdescribing a trajectory to be made, said electrophysiological signalsbeing preprocessed to obtain an observation tensor (X _(t)) at eachobservation time (t), changes in the observation tensor being modelledby a hidden Markov Model (HMM), said interface using a mixture of aplurality K of experts, each expert (E_(k)) being associated with ahidden state (k) of the HMM model and being defined by a multi-linearpredictive model, predictions of each expert being combined by means ofmixture coefficients (γ_(k) ^(t), k=1, . . . , K) to provide an estimate(Y _(t)) of a control tensor representing said command signals, wherein:said calibration is made during a plurality of calibration phases, atpredetermined instants of said trajectory, each calibration phase (u)corresponding to a plurality (ΔL) of successive observation times, andmaking use of a tensor X _(u) representing an observation tensor at saidplurality of successive times, a tensor Y _(u) representing a setcommand tensors at these same times and a matrix Z_(u) giving states ofthe HMM model at these same times, and in that each calibration phasecomprises: a step (a) in which tensors X _(u) Y _(u), are extractedusing the matrix Z_(u), observation tensors X _(u) ^(k) and commandtensors Y _(u) ^(k), k=1, . . . , K, relative to the different states(k) of the HMM model; a step (b) in which for each expert E_(k), tensors

and

are calculated corresponding to tensors X _(u) ^(k) and Y _(u) ^(k)respectively, after being centred and normalised, then a covariancetensor of tensor

and a cross-covariance tensors of tensors

and

being modified by adding covariance and cross-covariance tensors

and

respectively derived from a previous calibration phase, weighted by aforget factor, λ^(k); a step (c) in which the multi-linear predictiveexpert models are updated E_(k), k=1, . . . , K using a multivariateregression of partial least squares with exponential recursive weighting(REW-NPLS) regression, starting from covariance and cross-covariancetensors modified in the previous step; a step (d) in which the tensor

corresponding to the tensor X _(u), after being centred and normalised,the covariance tensor of tensor

and the cross-covariance tensor of

_(u) and Z_(u) are calculated, the covariance tensor of tensor

and cross-covariance tensor of

_(u) and Z_(u) being modified by adding to them the covariance tensor of

_(u−1) and the cross-covariance tensor of

_(u−1) et Z_(u−1), derived from the previous calibration phase, weightedby a forget factor λ; and a step (e) in which a second multi-linearpredictive model is updated giving a state vector ({circumflex over(z)}_(t)) of the HMM model as a function of the centred and normalisedinput tensor,

_(u), the components of the state vector providing possibilities thatthe interface is in each of the different states k=1, . . . , K at anobservation time, said update being made by an REW-NPLS regressionstarting from covariance and cross-covariance tensors modified in theprevious step, the mixture coefficients (γ_(k) ^(t), k=1, . . . , K) ofdifferent experts then being obtained using the second multi-linearpredictive model.
 2. The method of calibrating a direct neural interfaceaccording to claim 1, wherein in each calibration phase, a number oftransitions is determined for each pair of states observed during thiscalibration phase, and the transition matrix of the HMM model is updatedusing these numbers and the transition matrix of the HMM modelcalculated during the previous calibration phase, weighted by a forgetfactor.
 3. The method of calibrating a direct neural interface accordingto claim 1, wherein in step (b), a modified covariance tensor cov(

_(u) ^(k),

_(u) ^(k)) of the observation tensor for state k, centred andnormalised,

_(u) ^(k), is calculated by cov(

_(u) ^(k),

_(u) ^(k))=λ^(k) cov(

_(u−1) ^(k),

_(u−1) ^(k))+

_(u) ^(k)x₁

_(u) ^(k) in which x₁ is a tensor product for a first mode, and in thatthe modified cross-covariance tensor cov(

_(u) ^(k),

_(u) ^(k)) in which

_(u) ^(k) is a command tensor related to the centred and normalisedstate k is calculated by cov(

_(u) ^(k),

_(u) ^(k))=λ^(k) cov(

_(u−1) ^(k),

_(u−1) ^(k))+

_(u) ^(k)x₁

_(u) ^(k).
 4. The method of calibrating a direct neural interfaceaccording to claim 3, wherein in step (c), for each expert E_(k), k=1, .. . , K, an optimal rank, f_(opt) ^(k), in a latent variables space forthe REW-NPLS regression is determined as that which minimises an errorerr_(u) ^(k,f) ^(k) =λ^(k) err_(u−1) ^(k,f) ^(k) +∥Y _(u) ^(k)−Ŷ _(u−1)^(k,f) ^(k) ∥ in which err_(u) ^(k,f) ^(k) is a prediction error of thecommand tensor related to state k, by the expert E_(k) during thecalibration phase u for a rank f_(k)=1, . . . , F in which F is apredetermined maximum rank, and Ŷ _(u−1) ^(k,f) ^(k) is a prediction ofY _(u) ^(k) by the expert E_(k), obtained by the multi-linear predictivemodel of the same expert during the previous calibration phase, for thesame rank.
 5. The method of calibrating a direct neural interfaceaccording to claim 1, wherein in step (d), a modified covariance tensorcov(

_(u),

_(u)) of the observation tensor, centred and normalised,

_(u), is calculated by cov(

_(u),

_(u))=λ·cov(

_(u),

_(u))+

_(u)x₁

_(u) in which x₁ is a tensor product for a first mode, and in that thecross-covariance tensor cov(

_(u), Z_(u)) is modified by cov(

_(u), Z_(u))=λ·cov(

_(u), Z_(u))+

_(u)x₁Z_(u).
 6. The method of calibrating a direct neural interfaceaccording to claim 5, wherein in step (e), an optimal rank, f_(opt), ina latent variables space for the REW-NPLS regression is determined asthat which minimises an error err_(u) ^(f)=λ·err_(u−1) ^(f)+∥Y _(u)−Ŷ_(u−1) ^(f)∥ in which err_(u) ^(f) is a prediction error of the commandtensor during the calibration phase u for a rank f=1, . . . , F in whichF is a predetermined maximum rank, and Ŷ _(−u−1) ^(f) is a prediction ofY _(u) by a mixture of experts weighted by mixture coefficients obtainedby the second multi-linear predictive model during the previouscalibration phase, for the same rank.
 7. The method of calibrating adirect neural interface according to claim 1, wherein theelectrophysiological signals are electrocorticographic signals.